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, Abstract 

> 

We investigate entropy production for an expanding system of particles and resonances with isospin symmetry 
■ - in our case pions and p mesons - within the framework of relativistic kinetic theory. A cascade code to simulate 

the kinetic equations is developed and results for entropy production and particle spectra are presented. 

O 

1 Introduction 

' The search for the Quark-Gluon-Plasma (QGP) requires a detailed and quantitative understanding of the influence 
of known physics on the particle spectra and abundances from heavy ion collisions, such that possible systematic 
O ' deviations due to the new physics of the QGP can be recognized in the experimental data. 

When interpreting the pion momentum spectra JlJ-Q the question of entropy production during freeze-out occurred 
0| . If the contribution to entropy production after freeze-out due to resonance decays is known, the entropy of the 
^ : system before freeze-out can be calculated by measuring the multiplicities and momentum spectra of all particles after 
freeze-out and combining these data with an estimate of their initial spatial distribution (e.g. from a fireball model). 
It is then crucial to know how much entropy is produced in an expanding system during the transition (freeze-out) 
d from a highly interacting (equilibrium) state to a non-interacting one. Here we will give such an estimate for a simple 
system of pions and p mesons, using the framework of relativistic kinetic theory ||. This will enable us to estimate the 
reliability of models based on hydrodynamics in which freeze-out is usually implemented by cutting off the entropy- 
conserving expansion sharply at the so-called freeze-out surface, neglecting possible entropy production during the 
decoupling process or afterwards by the decay of unstable resonances. 

The paper is organized as follows: In Section || we will formulate via relativistic kinetic theory the problem of 
resonance formation and decay for a system of pions and p mesons, including Bose statistical effects. We will show 
that the collision term satisfies Boltzmann's H-theorem and discuss the formation and decay rates for p mesons. In 
Section 2 we describe a simulation of the kinetic equations of our system with the help of a cascade code, and in 
Section 4 we present our results, in particular on the production of entropy by resonance formation and decays. We 
will also discuss the influence of the interplay between p-formation and -decay on the time evolution of the particle 
spectra. In Section || we summarize our results. Some technical details are given in the Appendix. 
In the following we will use natural units, h = c = = 1, if not stated explicitly otherwise. 

2 Kinetic theory for the 7r-p-system 

As a starting point we consider the kinetic equation 

m^-f(x,p) = -T loss (x,p) f(x,p) +r gain (x,p) (l + f(x,p)), (1) 
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where r is the proper-time of a particle at point x with momentum p, and the right hand side of Eq. (|l]) is the collision 
term. The latter consists of a term for particle loss (with rate ri oss ) and one for particle gain (with rate r ga j n ). A 
vanishing collision term means vanishing entropy production (Liouville theorem). The left hand side of the kinetic 
equation can be rewritten as 

dr \ dr J dx» \ dr J dp^ dap dpt 1 

In the following we neglect external or mean-field forces by setting F 11 to zero. 

For 2 — > 2 processes Boltzmann's H-theorem can be proven from unitarity and particle number conservation alone 
H whereas in 2 ^ 1 reactions as in resonance formation and decay a different approach has to be taken. In || the 
bilateral normalization condition was exploited. We will here consider the isospin symmetry of the matrix elements 
to evaluate the gain and loss rates for the 7r-y0-system and thus prove the Boltzmann H-theorem explicitly. 

At low relative momentum the 7T7t interaction is dominated by the p resonance, with a small non-resonant back- 
ground (hJ. We are interested in the very last stage of a nuclear collision, after the pions and p mesons have "frozen 
out" , i.e. stopped interacting with other particles and with each other. After this freeze-out we expect the dominant 
physical process to be the decay of the unstable p mesons which we want to describe by kinetic theory. But in order 
for such a description to be consistent and the calculation of entropy production to be meaningful, our collision term 
should allow for detailed balance, which implies that we also have to include the inverse process tttt — > p (resonant 7T7T 
scattering), however small its rate should be. 

The interaction Hamiltonian for the 7r-p-system is of the form [[lO] 

H int = fpmT P M (<p x dptp) , (3) 
where the vectors refer to the isospin structure, and p M is a massive spin-1 vector field satisfying the gauge condition 

W = 0- (4) 

Due to energy momentum conservation and isospin symmetry only the processes p a «-> tt^tt 7 are allowed where 
a, /?, 7 € {0, +, — } are pairwise different. 

To write down the rates in terms of the distribution functions, we assume that the ensemble averages of multiple 
products of creation and annihilation operators arising from treating H ln t in first order perturbation theory can be 
factorized into products of ensemble averages of 1-particle operators, given by the following replacement (for p mesons 
a is a double index denoting both the spin and isospin component of the particle): 

{a a \k)a a ' \k')) = 5{k - k') f a (k) S aa ' (5) 
{a a {k)a a '\k')) = 6(k-k')(l + f a (k))5 aa ' (6) 
(a a (k)a a '(k')) = (a a Hk)a a 'Hk')) = 0. (7) 



In global thermal equilibrium this is justified by the thermal Wick theorem |jll|, and in 1 12 it was shown how to 
generalize the latter to local thermal equilibrium systems. For non-equilibrium ensembles the validity of this ansatz 
has been proven from first principles only for the </> 4 -theory (l3). We use it here because it leads to the classically 
expected result. 

Let us now introduce pion and p meson distribution functions f a = / 7rQ and F a = f Pa . In the case of the p, F a 
describes a single isospin and spin projection, but we assume spin saturation such that the distribution functions for 
the three spin projections are equal, and therefore suppress the spin index on F a . The loss rates are then given by 

rCM = (^f^Y^M I D P2Dp3(M( Pl ,p 3 )+M(p 2 ,p 3 ))Mp 2 )(l + F^(p 3 ))5(p 1 +p 2 -p 3 ), (8) 

0,1 



rf ° s (P3) = (2^) 4 I Dp 1 Dp 2 M(p 1 ,p 3 )(l + f- P ( Pl ))(l + f- ( (p 2 ))S(p 1 +p 2 -p 3 ), (9) 

Sp{S p -t- 1) ^ j 

while the gain rates are simply given by replacing f a ^> 1 + and F a «-> 1+F a . Here e Q /37 is the totally antisymmetric 
tensor, Dp is the Lorentz-invariant measure in momentum space 



where E = \fp 2 + m 2 , and the distribution functions are always taken on the mass-shell. 

The two contributions to the pion loss term in Eq. (g) arise from the two possibilities in (^) to let the derivative 
act on either one of the two scattering pions. 
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The momentum dependent part of the matrix element as obtained from (S) in first order perturbation theory is 
given by 



M(pi,p p ) = -pj + 



{PiPpF 1 (If, , x2 

-dr - = ^2(4 m±p P ) 



Pi-Po 



2 2 

P l m 



(11) 



••-p ••"p 

Here we have already summed over the three spin states of the p. To describe a single p spin state in a spin saturated 
system, Eqs. (^Tj) and (12) below must be multiplied by a factor l/s p (s p + 1) = 1/3, as in Eq. (^|). Since in Eq. (^) 
this factor is missing, this rate contains already the sum over all final p spin states . 

The dependence of the matrix element on the pion momentum pi arises from the derivative in the p7T7r-coupling 
([5]) , whereas the one on p 3 — p p comes from the p-polarization tensor |Tij which appears after summing over all three 
spin states. If the meson masses for different isospin states are taken to be equal, (ll|) can be simplified by using the 
momentum conserving (5-functions in (||j9[) to evaluate pi ± pj and exploiting the mass-shell conditions pf — mf : 



M(pi,p 3 ) = M(p 2 ,p 3 ) = 7 (r 



(2m*) 



(12) 



In this case the two contributions to the loss rate in (||) become equal. 
In kinetic theory the entropy current for bosons or fermions is given by 

Sp,(x) = / Dppf,, fi(x,p)]xifi(x,p) =F 0-±fi{x,p))\n(l±f i (x,p)) 

+fi(x,p) \nf t (x,p) T (1 ± Ji(x,p)) ln(l ± fi(x,p)) , 
with the upper (lower) sign for bosons (fermions) . The rate of entropy production can then be written as 

= -El d p [(^(^))in (j^-) + frrfap)) m (.Jl^ 



(13) 



(14) 



In our 7r-p-system all particles are bosons, i.e. the plus signs apply and there are no anti-particle contributions. The 
sum over i in (|l4| ) goes over the three charge states of the pions and p mesons and additionally over the three spin 
projections of the latter. After replacing p fJi d ,1 f(x,p) in (|l4) by the r.h.s. of the kinetic equation (|l|), the sum over 
the three p spin states cancels the factor l/s p (s p + 1) in Eq. (^), and one gets, after some renaming of integration 
variables and permutation of the isospin indices a, /3, 7, 



= (2^ 4 /p 2 7r ^ko/3 7 iy^Pl^2^3X(pi,P3)<5(pi+P2-p3) 

a/37 

x (l + f a (pi)) (l + fp(p 2 ))F^(p 3 ) (y - 1) lny , 



where 



f a (pi)ff}(p2)(l + F^(p 3 )) 



(15) 



(16) 



y (l + f a ( Pl ))(l + f p (p 2 ))F^(p 3 ) ■ 

y and M. are positive by definition. Since for y > the expression (y— 1) lny is positive semidefinite, we have d^S^ > 
always. This proves the Boltzmann H-theorem in our case and for similar systems with symmetric isospin structure. 
The collision integral ( fl5"| ) vanishes only if y = 1. 

Expressing (similarly to Eq. (|l3|)) the isospin current and the energy momentum tensor through the distribution 
functions fi(x 7 p) one can use the conservation laws for these quantities together with the condition for vanishing 
entropy production, y — 1, to derive the local equilibrium Bosc distribution 



Mx,Pi) 
1 + fi(x,Pi) 



exp 



(- \p( x )qi + u u (x)p"} /T{x) 



or 



fi(x,Pi) = 



,[n(x)qi+u u (x)pZ]/T(x) 



(17) 



where = 0, ±1 is the electric charge of species i. 

Important for a simulation of the kinetic equations will be the decay rate of the p meson. In the rest frame of a p 
with a given spin polarization it can be cast in the form 



loss 



( m p) 



2 

pTZTT 



l £ a/3 7 | 



167TS p (Sp + 1) 



d(a»0)(l + / /J (£))(l + / 7 (-jS)). 



(18) 



free decay rate rf ree 



The prefactor is the free decay rate rf roe , and the integral over the pion emission angles relative to the velocity of the 
p in the local heat bath gives the medium corrections due to Bose stimulated decay. The pion distribution functions 
under the integral are evaluated at \p\ = * m 2 — im 2 ,, the pion momentum in the rest frame of the decaying p meson. 
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In the limit of free decay in the vacuum the kinetic equation for p mesons takes the following simple form 



d 

Up 

This is solved by 



E p j t F a {p P ,t) = -T flcc F a ( Pp ,t) . (19) 



F a (p p ,t) = F tt (p p ,i ) e -( r ^/^)(*-M (20) 
The measured decay rate in the rest frame of the p meson (E p = m p ) is then given by 



f 2 (m 2 - 4m 2 ) 3 / 2 

rr ay K)-iF p 7 ■ ( 21 ) 

487T mi 



p 

Its experimental value at the maximum of the resonance at WL, = 769.9 ± 0.8 MeV is 151.2 ± 1.2 MeV. Together 
with the measured value of m n ± = 139.56995 ± 0.00035 MeV |1(| this determines the coupling constant as f p7T7I — 
6.049 ±0.027. 

For the simulation of the collision term we need a collision criterium which involves the 7T7r-scattering cross section 
tt(Pi) + 7r(p2) —* p{Pp)- Assuming a Breit-Wigner mass distribution for the p meson 

1 to r(m 2 ) f°° 

"^°» K- m y + ^ K ) ' L wim > m ' = 1 - <22) 

with a width given by the free decay rate r^^TO 2 ), Eq. ([H]), one gets, following standard textbook methods []l7|, |l8|, 

Mp (Mp ~ ^A) 2 ( , r , , \ l:m ..)_( jmc\ f , 

where 



a(j4; Pp ,l) = 4™ p (a, + I)/ 2 (1 _ ' / + ' _^ )3 + + »)) 



197.327 MeV 2 



?/),,, 



firi (23) 



/= 1fi f 7\,y Mp = ^, ^ = S L - ( 24 ) 

167TS p (s p + 1) TO p TO p 

Again we have summed over all three spin states of the p, and the factor (1 + F^(pi + P2)) in ( p3[ ) gives the Bose 
stimulation effect of the medium on p production. In an isospin symmetric system with F+ = F- = Fo it does not 
depend on the final isospin state 7, but only on the momentum p p — pi + pi of the p. 

The cross section and decay rate in vacuum and the mass distribution are shown in Figure |l| as functions of p p . 
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Figure 1: Free decay rate r p ccay (p; p ), free cross section a(p 2 ), and Breit-Wigner mass distribution w(jjl p ) = 2 H i p W(p, 2 ) 
for the p mesons as functions of their mass in units of to p , the mass at the maximum of the resonance. The position 
of the maxima of cr(p P ) and w(p p ) are (almost) identical, but at higher masses a has a lower weight than w(fj, p ). 

Finally we want to note that from a formal integration of the kinetic equations over the time parameter r it can be 
shown that a system with spherical symmetry in x and p remains spherically symmetric for all times (see Appendix 

S). 
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3 Simulation of the kinetic equations 

As initial conditions for the solution of the kinetic equations we start with a system of p's and 7r's in global thermody- 
namic equilibrium at vanishing chemical potential. The thermodynamic conditions will be chosen near the expected 
"freeze-out" point. The particles are distributed homogeneously and isotropically in a spherical box in coordinate 
space with an equilibrium momentum distribution (exp(E /T) — 1) . The masses of the p's are distributed with the 
Breit-Wigner distribution Eq. ( p2[ ) . For numerical reasons we will employ a cutoff for the p mass and for the momenta. 

We assume isospin symmetry in the initial state, and since this symmetry is conserved by the strong interactions, we 
can exploit it to reduce the number of coupled kinetic equations to be solved. We can sum the three equations for the 
pion distribution functions f a and the nine equations for the p meson distribution functions F as , where a = +, — , and 

s = +, — , denote isospin and spin projections, respectively. Setting /+ = /_ = fo = f and F ++ = = . . . = F. 

i.e. the distributions for all isospin and spin components equal to each other, we find the two coupled equations 

P p d p f(x,p) = -riJx,p)f(x, P )+ri ain (x,p)(l + f(x,p)), (25) 
p ti d t ,F{x,p,m p ) = -r^ os Jx,p,m p )F(x,p,m p )+r p gaiD (x,p,m p )(l + F(x,p,m p )) . (26) 

In the global coordinate system all distribution functions and rates are isotropic in momentum space and thus only 
functions of E. The rates are given by 

Tl ss (E) = 2^/ dm p w(m p ) JT / dE' f(E')(l + F(E + E\ m p )) , (27) 



r- ain (£) = 2^5/ dm p w(m p )^=== dE'(l + f(E'))F(E + E',m p ), (28) 



4tt J 2m „ " " p 'j&-ml J E - 



f2 A A pE+ 



KJE,m p ) = Jl / dE'(l + f(E'))(l + f(E-E')), (29) 

6 4?r JE 2 - ml Je- 



p 

2 



3 47r JE 2 - m 2 Je- 



with M. being the spin-summed matrix element (|l2]). The integration limits are 



E± = 

2m, 

E± = \{ E ^r\l ) * 

and w(m p ) = 2m p W(rn p ) is the mass distribution, Eq. (p2[), for the p mesons. The distribution functions are normalized 
according to 

N„ = N„+ + N„- + N„c = 3 J ^iff(r, P, t) 

= 3^ J r 2 dr J dE E y/ E 2 - mjf(r, E, t) , (32) 



i d^v d^p i 

N p = N p + + N p - +N p o = 9 J 3 I dm p w(m p )F(r,p,m p ,t) 



9- [ r 2 dr [ dm p w{m p ) [ dEE^I E 2 - m 2 F(r, E, m p , t) , (33) 
J J2m, J m „ V ^ 



with N p + N n /2 = N being constant during the time evolution. 

In the entropy production rate, Eq. (15), all 6 terms of the sum become equal, and it simplifies to 

0% = 6(27r)V£™- / dm p w(m p )M{m p ) J Dp 1 Dp 2 Dp 3 S(j) 1 +p 2 -p 3 ) 

f(pi)f(P2)(l + F(p 3 ,m p )) \ 



x In 



1 + f(pi))0- + f(p2))F(j>3,m p ) J 
(/(Pi)/(P2)(1 + F(p 3 ,m p )) - (1 + /(pi))(l + f(p2))F(p 3 ,m p )) . (34) 
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With the help of the four-dimensional <5-function the p 3 and m p integrations can be evaluated. This leads to 

M{m p ) 



d»S u 



6tt/, 



piXTT 



x In 



Dp 1 Dp 2 w(m p ) ■ 

trip 

' f(pi)f(p 2 )(l + F(p 3 ,m p )) 



x [l + f(j?i))(l + f(p2))F(p 3 ,m p ) J 
x (/(Pi)/(f> 2 )(l + F(p 3 , m p )) - (1 + f(pi)) (1 + f(p2))F( P3 , m p )) , 
with p3 = pi 



(35) 



P2, 



= V / (Fl+F 2 ) 2 -(p 1 +p 2 ) 2 . 



To simulate the time evolution of the particles we wrote a cascade code in which the distribution functions / and F 
are represented by isospin symmetric test particles. In terms of the initial radius R of the system the conserved particle 
number N is given by N = g ■ R 3 , where g involves the initial ir and p densities and is a function of temperature, 
see Eqs. ( |32| , (33|) . The initial radius of the fireball was fixed to be 5 fm; the corresponding particle numbers for 
temperatures T = 100, 150 and 200 MeV are given in Tabic g. 



Table 1: Particle numbers, density g, and the conserved particle number N for a system of p mesons and pions in 
equilibrium at T = 100, 150, 200 MeV, for a fireball radius of R = 5 fm. 



T [MeV] 






e [fm- 3 ] 


N = + N p 


100 


16.8 


0.6 


0.07 


9 


150 


61.8 


11.1 


0.34 


42 


200 


172.0 


70.0 


1.25 


156 



For the evaluation of the integrals over the p mass distribution a cutoff at 2.5 m p = 1920 MeV was introduced, 
and the momentum integrals were cut off at p m ax = 3.5m p = 2689 MeV. The mass cutoff was large enough to have no 
influence on our results in spite of the fact that for this cutoff the tail of the Breit-Wigner distribution still contains 
about 10% of the mass spectrum. The reason for this is the exponential form of the bosonic energy distribution in 
Eq. ( |33| ) which suppresses the formation of high-mass p's and results in a deviation of the mass distribution for thermal 
p's from the mass distribution of p's at rest (see Fig. |J). Thereby the peak of the mass spectrum and its weight is 
shifted towards lower masses as the temperature decreases. A second peak close to the threshold of p production from 
two pions appears for temperatures below T « 75 MeV. It dominates the spectrum for low temperatures, but since 
we take the 7r's and p's to be in local equilibrium the number of p mesons becomes very small for these temperatures 
and vanishes for T = 0. 

The code determines simultaneously the space and momentum coordinates of all test particles as functions of a 
global time parameter on a hyper-surface in phase space. (For an introduction to (nuclear) cascade codes see ]Ts||,po|; 
the problem of a simultaneous relativistic time parameter is discussed for instance in p^| .) Since no mean field 
effects are taken into account, the particles propagate freely on straight lines between p formation and decay: 

x (r + At) = x(t) + (3 At, P=^, 

p(t + At) = p(r). (36) 

r is the global time parameter. 

The formation of a p meson in the scattering of two pions is implemented via the collision criterion 

Tib 2 <2a(m 2 p ,E p ). (37) 

The factor 2 arises from the same factor in front of the matrix element in Eqs. ( p7|,p8| ). <j{m 2 pl E p ) is the medium- 
modified cross section (p3|), and b 2 is the relativistically invariant impact parameter J22| 

^-( Al -^p) 2 = - (Al) = + (f^Q!. (38) 

Here P = p\ + p 2 is the momentum of the pion pair, and P 2 — m 2 is the square of their invariant mass. In the 
center-of-mass system of the two pions the impact parameter becomes b 2 = Ax 2 m - If b satisfies the collision criterion 
Eq. (|37|) the collision of two pions takes place in the global system at the point where their distance 

- (Ax) 2 = (Ax) 2 = (Aa? (r ) + A/3 • r) 2 , A/3 = || - g (39) 
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500 1000 1500 

Mass m p [MeV] 

Figure 2: Mass distributions for p mesons in equilibrium with pions for different temperatures. The temperatures from 
top to bottom are T = 200, 175, 150, 125, 100, 75, and 50 MeV. Also shown for comparison is a pure Breit-Wigner 
mass distribution (dotted line), normalized at its peak to the T = 200 MeV curve. 



is minimal. This happens at the collision time t c which is given by 

= Ax(r )-A/3 
Tc (A/3) 2 ' 



(40) 



and which is measured relative to the last grid point to of the global time parameter at which all trajectories were 
calculated. 

For the decay of a p meson we assume isotropic decay in its rest frame. To decide if a decay in the rest frame of 



the p takes place we use the mass-dependent decay rate Eq. (18) which involves the pion distribution function in the 



rest frame of the decaying p meson. We evaluate the latter by counting pions in phase space cells around the pion 
momentum p, \p\ = \Jnif, — 4m 2 . To avoid the transformation of all pion coordinates from the global system into the 
rest frame of the decaying p meson, we transform instead the phase space volume into the global system and do the 
counting there. 

The phase space contributions to the cross section and decay rate are averaged over a large number of parallel 
ensembles. The latter influence each other only via the evaluation of the Bose factors as an ensemble average. To 
obtain from the discrete locations of the test particles smooth, but positive definite phase space density distributions, 
we smear out the contribution from a particle at some point {x,p) over phase space with a Gaussian weight function. 
According to the uncertainty principle the smallest possible phase space cell has the volume A 3 xA 3 p = h 3 , thus the 
widths of the weight functions in x and p have to fulfill AxAp = 2irh. 

In evaluating the collision term, including the Bose enhancement factors, with probabilistic methods, we use a 
trick similar to the one described in |2j| (see also |24|, |25|). Instead of computing for each possible p decay or creation 
the exact contributions from the Bose factors to the decay rate and cross section, we first determine the decay or 
formation probability using an upper bound for the Bose factors (and thus the transition rate) . Only if the collision 
would actually happen under these conditions, we proceed to check the collision criterium with a more realistic estimate 
of the Bose enhancement factors. This saves an appreciable amount of computer time by eliminating without effort a 
large number of "unsuccessful collision attempts" . The upper bound on the Bose factors is refined in two successive 
steps: In a first step we obtain a global (and time-independent) upper limit on the Bose factors from the initial phase- 
space density of the expanding system. If by using this global limit we find that p formation or decay are possible, 
we calculate a more accurate dynamical bound from the actual density distributions of the previous time step. Only 
if this second upper bound also allows p formation and decay do we actually calculate the distribution functions as 
described above. 

In our simulations the number of particles in each system, N = N 7T /2-\- N pi times the number of parallel ensembles 
was fixed to w 300000. We followed the system for 150 time steps of At = 0.1 fm/c. 
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4 Results 



4.1 Decay and re-scattering rates 

In general the decay of the p mesons is governed by the decay law 

N p {t + At) = -R 3 / dm p w(m p ) / dp p 2 F(E, m p , T )e- r >°- (£ ' m " )AT/B , (41) 

n J2m„ JO 

where E 2 = p 2 + m 2 and r^ ss is the loss rate for p mesons from Eq. (^). If we neglect the influence from Bose 
stimulation by surrounding pions we can replace rf oss /_B in the exponent by r^ ecay {m p )/j from Eq. ( pi] ) where 7 is 
the Lorentz factor between the global system and the p rest frame. 

For the decay of p mesons at rest (7 = 1) into the vacuum one then obtains 

N p (t + At) = J dm p N p (m p , t) exp(-r^ cay (m p ) Ar/ 7 ) , (42) 

where N p (m p , t) is the number of p mesons with mass m p at time t. Figure |^ shows a numerical simulation of this 
decay law with the help of our kinetic code. There is good agreement between the simulated result and the one 
calculated directly from Eq. (|^j). The striking deviation of both these results from the straight line describing the 
exponential decay of p mesons with mass m p = 770 MeV and a constant width Yf^ = 151 MeV arises from the mass 
dependent decay rate (|2l]): The heavier p mesons from the initial Breit-Wigner mass distribution decay first (due to 
their larger than average width the slope of the decay curve is at small t actually slightly steeper than the reference 
line), leaving at later times only the more long lived lighter p mesons from the lower end of the mass distribution. 
Therefore the decay curve levels off at large t. 




5 10 15 

Time x [fm] 

Figure 3: The number of p mesons as a function of time for free decay at rest simulated with 1000 parallel ensembles 
each containing initially 300 p mesons (solid) and calculated according to Eq. (ji|) (dashed) . For comparison we show 
the free decay for a constant decay width r^ ecay = 151 MeV (dotted). 

Fig. |J shows the p meson decay for initially thermalized systems with different initial temperatures when we neglect 
the back-reaction tttt — > p. One clearly sees a deviation from the just discussed decay of p mesons at rest, shown again 
by the dotted curve. A comparison between the dashed and dash-dotted curves in Fig. [| shows that this effect is not 
due to medium effects via Bose stimulation from the pions in the system. It is, however, rather accurately reproduced 
by an analytical calculation (solid lines) which uses Eq. ( f4l| ) with a thermal distribution function and again neglects 
Bose stimulation: 



NJt + At) = -W 



dm p w(m p ) 

2m„ ' JO 



dp p 



l/T 



- 1 



-rf"K)Ar/ 7 ( P ,m,,) 



(43) 



For all three initial temperatures the curves in Fig 
since the Lorentz-dilatation factor j(p, m p ) = Jp 1 



have approximately the same form. This was not expected, 



m 2 Jm p in the exponent of Eq. (H3J) is weighted quite differently 
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5 10 15 

Time x [fm] 

Figure 4: p decay for systems without p formation for different initial temperatures. Dashed and dashed-dotted: 
simulation with and without phase space contributions to the decay rate. Solid: calculated time evolution for the decay 
of p mesons with initial equilibrium distribution. From top to bottom: T = 200, 150 and 100 MeV. For comparison: 
center of mass decay (dashed-dotted) as calculated with Eq. and fitted to the initial number of p mesons for T — 
200 MeV. 



at different temperatures, which should lead to an effectively larger decay rate at lower temperature (with smaller 
average (7)). However, at lower initial temperatures also the mass spectrum of the p mesons is shifted towards lower 
masses (see Fig. |^). Since the p width decreases with decreasing mass m p , the net effect is a somewhat lower effective 
p meson decay rate at lower temperature. 

The observed weakness of medium effects from Bose stimulation shows that the chosen initial conditions select 
systems which are already quite dilute in phase space. Bose statistical corrections would play a larger role if the 
systems were not allowed to expand (2^, [24|, [25) . 

In Fig. H we study the additional effects from the back-reaction irir — ► p. From a formal point of view they must 
be included in order to preserve detailed balance and the Boltzmann H-theorem. In practice one sees that the number 
of resonances now decreases much more slowly (solid and dotted curves in Figs. ||a-c). Again, the inclusion of Bose 
enhancement factors in decay rates and scattering cross sections has no visible effect. 

Another illustration of the important practical role of p formation by pion re-scattering is shown in Figs. ^ where 
the total rate for the change of the number of p mesons is split into its loss and gain contributions. Since the system is 
initialized in thermodynamic equilibrium, the two contributions should initially balance each other. That in Fig. (0) 
they don't quite do that is due to the fact that at r = the walls around the initial sphere are removed and thus p 
mesons from the surface can now escape without re-scattering. One sees that initially the back reaction rate is still a 
considerable fraction of the decay rate, and that only for large times (t > 5 — 10 fm/c) the gain of p resonances by pion 
re-scattering can be neglected compared to the losses by decay. Furthermore, even if the absolute numbers for the decay 
and formation rates decrease drastically with decreasing temperature due to the exponentially decreasing occupancy 
numbers, the time dependence of the ratio between these two rates shows no strong temperature dependence. 

In order to better understand these features we give in Table | the pion mean free path Z£. 0G for different temper- 
atures, as extracted from the mean thermal velocity (v)t and the effective collision time: lf roe = {v)tTco\\- The mean 
thermal velocity of the pions at the beginning of the time evolution is given by 



(v)t 



^(p/E)[e E / T -ir^dp 

J™[ e E/T_ l]-l p 2 dp ■ 



(44) 



The pion (inelastic) collision time can be calculated from p meson gain rate, using the conservation of N = N p + 
in the process tt + ir — > p. Defining t* oVi = —N^/N^ = +N 7r / (2N P ) we get 



A^(t) 



2(A7V g " ain (r)/Ar) ' 



(45) 



with AN^ in (r) taken from Figs. | at r = 0. 
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Figure 5: p decay from our simulation for systems with different initial temperatures. Curves from top to bottom: 
p decay plus back reaction tttt — > p, with (solid) and without (dotted) Bose enhancement factors; p decay without 
back reaction, with (dashed) and without (dash-dotted) Bose enhancement factors; free decay rate (—•■ — ■■—) for 
comparison. 




5 10 15 5 10 15 5 10 15 
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Figure 6: Gain and loss rates of p mesons as a function of time for different initial temperatures. 



The system with T = 100 MeV has a pion mean free path which is much larger than the initial radius which in our 
simulations was chosen as R = 5 fm. If we assume the usual definition of "freeze out" which says, that particles will 
no longer hit each other if their mean free path becomes larger than the radius of the system, then "freeze out" will 
occur immediately after the beginning of the time evolution. In terms of absolute numbers for p formation (Fig. ^|c) 
this is fulfilled to good accuracy in our simulation. On the other hand the system with T = 200 MeV has a pion mean 
free path which is smaller than the radius of the system, so that, at least at the beginning of the time evolution, it 
will stay close to equilibrium. Nevertheless the shape of the curves for p decay in Fig. [| and for the differential gain 
and loss of p mesons in Fig. ^ are quite similar for the different temperatures, and show that the interaction itself 
does not change its character, only the number of particles involved gets smaller from T = 200 MeV to T = 100 MeV. 
Therefore "freeze out" is a quantitative statement about the number of possible interactions and not a qualitative 
statement about the form of the interaction. 
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Table 2: Collision time r^ oll , average thermal velocity (v)t, and mean free path lf Yac for the pions at the beginning of 
the kinetic evolution for different initial temperatures. The initial radius was 5 fm. 



T [MeV] 




(v) T 


Ifree [fm] 


100 


50.0 


0.811 


40.5 


150 


7.6 


0.867 


6.6 


200 


3.4 


0.899 


3.0 



4.2 Time evolution of the p mass spectrum 

The mass spectrum of the p mesons is changed both by decay and formation processes, see Fig. [?]. For the decay of p 
mesons at rest and for the decay of p mesons with an initial thermal momentum distribution we can calculate the time 
evolution for the mass spectrum from Eqs. ([|2]) and respectively. The results are shown as the smooth curves in 
Figs. 0b, c. The temporal change of the mass spectrum arises mostly from the mass dependent decay rate. We observe 
that the spectra for large times become qualitatively similar in shape to the initial spectra at low temperatures (see 
Fig. ||) and have their peaks shifted towards the p formation threshold at 2m 1T - Comparing Figs. and 0a we see 
that the effect of p formation by pion re-scattering is a reduction of the effective decay rate for the p mesons. If in 
the two different time evolutions we select points with an equal total number of p mesons, the shape of the p mass 
spectrum is, however, nearly identical. 




500 



1000 



1500 



500 
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1000 



1500 



1000 1500 

Mass m p [MeV] 

Figure 7: Mass spectrum of p mesons as a function of time from the kinetic simulation, for an initially thermalized 
system with temperature T — 200 MeV, in time steps of 3 fm/c. The Breit-Wigner distribution, normalized at its 
peak (vertical dashed line) to the initial thermally smeared mass distribution, is shown as a dotted line for comparison, 
a: total system including decays and pion re-scattering, b: p decays only, with the smooth solid curves indicating the 
result from Eq. (f43|). c: decay of p mesons at rest, with the smooth curves indicating the result from Eq. (E2J) . 



4.3 Single particle spectra 

The pion spectrum from decaying p mesons at rest is fully determined by the initial p mass distribution and the 
two-body decay kinematics. This is shown in Fig. || The maximum of the pion energy spectrum at E v = m p /2 = 385 
MeV results from the maximum of the resonance at rn p = 770 MeV, 

In the absence of pion re-scattering, in our model the Lorentz-invariant pion momentum spectrum consists of a 
superposition of the initial thermal pion spectrum and the spectrum of decay pions from thermally distributed p 
mesons M . This is shown as the solid line in Fig. where the contribution of pions from direct p decay (dotted line) 
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is given by 



(27T) 



i dN I _ 
d 3 p 



9- 



1 



Pi/p 2 + to 2 
™2 r 



dm to w(m) 



E dE 



4m2 J e? E 

E- 



(46) 



E" 



2mi 



\/p 2 + ml ± — v/rf 

TO 



Aral 



Once p formation by pion re-scattering is taken into account, one might expect that this changes both the direct pion 
spectrum (because pion loss by re-scattering into the p channel is not equally distributed in momentum space) and the 
decay spectrum from p decays (because the p's from re-scattering pions need no longer have a thermal distribution) 
after the system is allowed to expand. (Of course, initially the decay and re-scattering processes don't change either 
spectrum because the system is in thermodynamic equilibrium.) Figs. || show that the continuous nature of freeze-out 
in our kinetic simulation has only very weak effects on the final energy spectra: even for the rather dense system 
with initial temperature T — 200 MeV the somewhat steeper final pion spectrum (compared to the initial thermal 
one) is very accurately described by a thermal direct pion component superimposed by decay pions from thermally 
distributed p's. (We have checked that our simulation, with re-scattering shut off such that all p mesons decay directly 
with their initial thermal distribution, agrees with the resonance decay calculations of Eq.(f46|) .) No additional cooling 
of the pions by the re-scattering processes is observed, if it is there it is accurately compensated by the developing 
radial collective expansion. We conclude that the calculation of the single particle spectra according to the methods 
used in Refs. ||, |], [|, ^| (which do not include any kinetic evolution after the so-called freeze-out point) provides a 
quantitatively accurate approximation. 



4.4 Entropy production during expansion 

Before attempting to calculate the entropy production in our system it is very important to gain a rough qualita- 
tive picture of the phase space distributions of the particles. At the beginning the distributions in momentum and 
coordinate space are totally uncorrelated, due to the assumption of (global) thermal equilibrium. After allowing the 
system to expand, strong correlation between the momenta and coordinates develop. After a short interactive phase, 
the system essentially evolves by free-streaming, and after some time the momentum distribution at each point in 
coordinate space becomes essentially a ^-function in momentum space 8(r — (p/E)t). This very uneven distribution 
in phase-space leads to severe problems if one tries to calculate the entropy by a naive phase-space integration of 
/ In / without violating numerically Liouville's theorem that the entropy in system of free-streaming particles remains 
constant. 

We approached this problem by a particular binning in phase space adapted to a system of free-streaming particles. 
We evaluated the total entropy and the entropy production rate on phase space cells of fixed size Ax • Ap = 2irh, 
whose actual shape, however, was adapted to the free expansion by projecting the particle coordinates for each time 
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Figure 9: Pion momentum distribution as a function of pion energy in the global reference frame at the end of the time 
evolution. The circles indicate the results from the numerical solution, which includes both p decay and formation 
by pion re-scattering. The error bars indicate the statistical error from the sampling over parallel ensembles. For 
comparison: equilibrium (Bose) pion distribution at the beginning of the time-evolution (dashed line) , and (as solid 
curve) a superposition of these initial pions with pions from p decays (dotted line) , assuming that all initially present 
p mesons decay directly, without pion re-scattering. 

step to the point in space where the particles would have started their trajectory if no interactions and decays had 
taken place. We then evaluated the particle densities within cubic boxes in the new coordinates, thereby obtaining an 
equivalent "initial" phase-space distribution for free streaming particles. 

This method is applicable as long as only a small additional volume in phase space is occupied by particle inter- 
actions during time evolution. This is the case for systems which freeze out shortly after the beginning of the time 
evolution. If the system spread out more into the initially available free phase space, the phase-space densities could 
become too low for a reliable simulation of the kinetic equations by our test particle method. 




Time x [fm] Time x [fm] Time x [fm] 

Figure 10: The total entropy as well as the individual contributions from pions and p mesons as a function of time. 
Dashed curves: p decay only, no pion re-scattering. Solid curves: Full numerical simulation including re-scattering 
processes. 
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In Fig. [10| we show the calculation of entropy as a function of global time r. It is calculated from Eq. (|13|), by 
integrating the entropy current over a hyper-surface of constant global time: 

S( T ) = J S»{x)dEfr ) (x) = J S°(x,T)d 3 x. (47) 

Eq. ( |l3|) allows to separate the contributions from pions and p's. The numerical simulation was stopped after t = 15 
fm/c when nearly all p mesons had decayed and the entropy had saturated. Table [| gives the initial equilibrium values 
for the entropy, Table ^ the final values and the relative increase during the time evolution, both for p decays without 
re-scattering and including re-scattering. 



Table 3: Initial equilibrium values for the entropy for different initial temperatures. 



T[McV] 


s P 




Stot 


200 


469.8 


738.4 


1208.2 


150 


80.3 


301.1 


381.5 


100 


4.3 


90.1 


94.4 



Table 4: Final entropy and relative entropy increase for three different initial conditions. 



a) p decay only 


T[MeV] 


Sp 


Sir 


Stot 


AS/Stot in % 


200 


4.2 


1246.8 


1251.1 


4 


150 


0.8 


391.2 


392.0 


3 


100 


0.0 


95.1 


95.1 


1 


b) p decay plus pion re-scattering 


T[MeV] 


Sp 


Sir 


Stot 


AS/Stot in % 


200 


17.9 


1369.0 


1386.9 


15 


150 


1.9 


410.9 


412.7 


8 


100 


0.1 


96.1 


96.2 


2 



One sees that p decays alone produce very little additional entropy. Even for an initial temperature of T = 200 
MeV, where at the beginning nearly 30% of all particles are p mesons (see Table Q), the total increase in entropy is 
only 4% of the initial value. Including the re-scattering processes leads to a somewhat larger entropy increase, by up 
to 15% at T = 200 MeV. One sees from the numbers that in general the re-scattering processes are much more efficient 
in producing entropy than the decays. But even at the highest temperatures investigated in this work the entropy 
increase is modest; for systems which are initialized near freeze-out (T < 150 MeV, see Table ||) it is well below 10%. 
This implies that, as far as calculating the entropy balance from the final particle distributions is concerned 0, the 
fact that freeze-out is not a sudden process can be safely ignored. 



5 Conclusions 

We presented a cascade simulation with thermal equilibrium initial conditions for the kinetic evolution of the ir-p- 
system. This work was motivated by an earlier investigation of the influence of resonance decays on the entropy 
balance in relativistic heavy-ion collisions which did not take into account that particle freeze-out is a continuous 
rather than a sudden process. We have answered the question to what extent additional entropy can be created during 
the final kinetic evolution between the time when local thermal equilibrium is first broken and the time when finally 
all resonances have decayed. We found that this additional amount of entropy is very small, less than about 10% of 
the entropy already present at the beginning of decoupling. We also found that this additional entropy production is 
dominated by the resonant re-scattering between the pions in the system rather than by the resonance decay itself. 
The entropy saturates once the system is so dilute that re-scattering events become unlikely. 

Our simulations took into account the Bose statistics of the mesons, and we checked that our collision term satisfies 
Boltzmann's H-theorem and eventually leads to Bose equilibrium distributions if the matter is not allowed to expand. 
We realized, however, that in the calculation of the collision rates medium corrections by stimulated emission were 
very small. Consequently the influence of Bose statistics on the kinetics and dynamics of the expanding system turned 
out to be negligible. 
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We also found that the final pion single particle spectra resulting from the kinetic simulation can be very well 
approximated by the initial thermal pion distribution plus a contribution of pions from the decay of p mesons with a 
thermal momentum distribution with the same initial temperature. In other words, the modifications of the spectral 
shape from the final kinetic evolution of the 7r-p-system, in particular the effects of pion re-scattering into p-mesons, 
are minor. Thus the usual method pi ^| of calculating the resonance decay contributions to the pion spectrum in 
the "sudden approximation" , i.e. without considering the detailed kinetic balance between decays and re-scattering, is 
a good approximation and quantitatively quite accurate. This is certainly good news for practitioner because it means 
that for a thermo- and hydrodynamic analysis of measured hadron spectra the large numerical efforts of a detailed 
kinetic treatment of the freeze-out stage are not necessary. 

In the course of our investigations we also discovered a number of other interesting aspects of our system. It is well 
known that in a thermalized system the p mass spectrum is shifted towards lower masses compared to the vacuum 
case, due to the weighting with the exponential thermal Bose distribution (see Eq. (^3|)). However, during the time 
evolution we observed a further downward shift in the mass peak due to the mass dependent decay rate. This effect 
simulates very efficiently an increasingly cooler environment for the p mesons as time goes on. 
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initial stages of this work. U.H. wishes to thank B. Muller for stimulating discussions and the Physics Department at 
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A Appendix 



For a system with spherically symmetric initial conditions in p and x we show by formal integration over the global 
time parameter r that the kinetic equations conserve the symmetry. 
From the kinetic equations we obtain formally at time r + At 



F a (r + At, m p ) = F a (r, m p 



(2tt)< 



f 2 

J P7T7T 



m p s p (s p + 1) 



-r ^2 £a Pt / D PP D Pi s(pp + 

P7 Pa )M 



F Q (T,m p )(l + / /3 (r))(l + / 7 (T)) - (l + F a (T,m p ))fp(r)f 7 (r) 



At. 



(48) 



f 2 r 

/a(r + Ar)=/ a (T) + 2(27r) 4 ^^e Q/37 / D V(i D Pl S(p a + Pp - Pl )M 



(1 + Ur)) (1 + Mt))F 7 (t, m p ) - / q (t)/ /3 (t)(1 + F 7 (r, m p )) 



At, 



(49) 



where /(t) is short for /(x(t),p(t)) and 



x(t + At) = x(t) + ?^-At 

E 

F(t) 

p(t + At) = p(t) + — ^At =p(r). 

E 



(50) 
(51) 



With the above mentioned spherically symmetric initial conditions in space it follows immediately that / remains 
a function of only |af| for all times. To prove the same for the momentum dependence we have to show that the 
momentum integrations yield a result which is independent of the direction of p a . But the p a -dependent distribution 
functions F a , 1 + F a and f a , 1 + f a as well as the momentum independent matrix element ( |l^ ) can be pulled out of 
the integrals, and the the remaining integrations are completely independent of p a . Thus the spherical symmetry of 
the momentum distribution is also preserved. 
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